#include "OptData.h"
#include "LBFGS.h"
#include "LBFGSBopt.h"
#include "ConjuGradient.h"
#include "GeneralSetup.h"
#include "anyoption.h"
#include "MDopt.h"
#include "boost/date_time/posix_time/posix_time.hpp"

using namespace std;
using namespace boost::posix_time;

int main(int argc, char* argv[]) {
	AnyOption opt;
	opt.addUsage("");
	opt.addUsage("Usage: ");
	opt.addUsage("--setupFile                       setupFile file");
	//opt.addUsage("--lowerBound                      lowerBound of weights");
	opt.setOption("setupFile");
	//opt.setOption("lowerBound");
	opt.processCommandArgs(argc, argv);
	string setupFile = "";
	//double LBound=0;
	if (opt.getValue("setupFile") != NULL) {
		setupFile = opt.getValue("setupFile");
	} else {
		opt.printUsage();
		return 0;
	}
	//if (opt.getValue("lowerBound") != NULL) {
	//	LBound = atof(opt.getValue("lowerBound"));
	//}
	//cout<<"LBound = "<<LBound<<endl;

    GeneralSetup setup;
	setup.parseSetupFile(setupFile);

	OptData optData;
    optData.inputDoseSize(setup.dosegridFile);   // input the size of dose matrix

	for(int indexField=0;indexField<setup.fieldFiles.size();indexField++){
		optData.inputField(indexField,setup.fieldFiles[indexField].dijFile,setup.fieldFiles[indexField].weightFile);
	}

	for(int indexStruct=0;indexStruct<setup.structureFiles.size();indexStruct++){
		optData.inputStructure(setup.structureFiles[indexStruct].name,setup.structureFiles[indexStruct].structureFile,setup.structureFiles[indexStruct].eud);
	}
	
    vector<SConstraint> vConstrains; // add constrains 
	vConstrains = setup.objectives;

	vector<double> objFunValue;

	ptime start = boost::posix_time::microsec_clock::local_time();

    //optimization
    vector<double> resultVWeight;
    
	LBFGSopt *lbfgs = LBFGSopt::Instance(); // use LBFGS to solve the problem
	// traditional optimization  objective function is (DijWij^2-Dp)^2
    //lbfgs->optimization(&resultVWeight, &optData, &vConstrains, &objFunValue,setup.maxIter); 

	// traditional optimization  objective function is (DijWij^2-Dp)^2
	// program will check weights of beamlets to guarantee the initial value of weights used in algorithm is not equal to 0   					
	lbfgs->optimization_modify_1(&resultVWeight, &optData, &vConstrains, &objFunValue,setup.maxIter); 
																				
	//modified objective function, it is still (DijWij-Dp)^2
	//during the algorithm, program will check weights of beamlets not smaller than 0. If it is smaller than 0, it will be set to 0.   
	//lbfgs->optimization_modify_2(&resultVWeight, &optData, &vConstrains, &objFunValue,setup.maxIter); // traditional optimization  objective function is (DijWij^2-Dp)^2

	//ConjuGradientOpt *cg = ConjuGradientOpt::Instance(); // use Conjugate Gradient to solve the problem
	//cg->optimization(&resultVWeight, &optData, &vConstrains, &objFunValue,setup.maxIter); // optimization

	//MDopt mdOpt;
	//mdOpt.optimization(&resultVWeight,&optData,&vConstrains,&objFunValue,setup.maxIter);

	//LBFGSBopt *lbfgsb = LBFGSBopt::Instance();
	//lbfgsb->optimization(&resultVWeight, &optData, &vConstrains, &objFunValue,setup.maxIter,LBound); 


	ptime stop = boost::posix_time::microsec_clock::local_time();
	time_duration diff = stop - start;
	long milliseconds = diff.total_milliseconds();
	std::cout << milliseconds << std::endl;
	objFunValue.push_back(milliseconds);
	
    //output result data
    optData.outputWeight(&resultVWeight, setup.outputDir);
    vector<double> resultDose;
    resultDose = optData.calculateDose(&resultVWeight);
    optData.outputDose(&resultDose,setup.outputDir);
	optData.outputObjValue(&objFunValue,setup.outputDir);
	optData.outputDVH(&resultDose, setup.outputDir);
	
    return 0;
}